NATRE microstructure dataset#

Study NATRE dataset, reproduce Ferrari & Polzin (2005)

%load_ext watermark
%matplotlib inline

import glob

import cf_xarray as cfxr
import dcpy
import distributed
import gsw_xarray
import hvplot.xarray
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import tqdm
import xfilter
from IPython.display import Image

import eddydiff as ed
import xarray as xr

xr.set_options(keep_attrs=True)

plt.rcParams["figure.dpi"] = 140
plt.rcParams["savefig.dpi"] = 200
plt.style.use("ggplot")

%watermark -iv

/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/statsmodels/compat/pandas.py:65: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  from pandas import Int64Index as NumericIndex

numpy      : 1.22.3
tqdm       : 4.64.0
distributed: 2022.4.0
hvplot     : 0.7.3
cf_xarray  : 0.7.1.dev16+g7305fd8
gsw_xarray : 0.3.0
xarray     : 2022.6.0
matplotlib : 3.5.1
eddydiff   : 0.1
dcpy       : 0.2.1.dev11+gebbbf87.d20220428
xfilter    : 0.2.1.dev6+g5fdc006

if "client" in locals():
    client.cluster.close()
    client.close()
client = distributed.Client(n_workers=2, threads_per_worker=2)
client

Client

Client-8948807e-35dd-11ed-95d3-3af9d394f1c6

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Cluster Info

LocalCluster

f257c023

Dashboard: http://127.0.0.1:8787/status	Workers: 2
Total threads: 4	Total memory: 16.00 GiB
Status: running	Using processes: True

Scheduler Info

Scheduler

Scheduler-35410ff6-70d8-4bb7-85bd-0a821e8abd2a

Comm: tcp://127.0.0.1:54020	Workers: 2
Dashboard: http://127.0.0.1:8787/status	Total threads: 4
Started: Just now	Total memory: 16.00 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:54032	Total threads: 2
Dashboard: http://127.0.0.1:54034/status	Memory: 8.00 GiB
Nanny: tcp://127.0.0.1:54024
Local directory: /Users/dcherian/work/eddydiff/notebooks/dask-worker-space/worker-c83c9tud

Worker: 1

Comm: tcp://127.0.0.1:54031	Total threads: 2
Dashboard: http://127.0.0.1:54033/status	Memory: 8.00 GiB
Nanny: tcp://127.0.0.1:54023
Local directory: /Users/dcherian/work/eddydiff/notebooks/dask-worker-space/worker-5g_85fsd

Image("../images/natre-large-scale.png")

Read dataset#

created using eddydiff.natre.combine_natre_files()

natre = ed.natre.read_natre(load=True)

13 40
0  values == 0
28090

/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/numba/np/ufunc/gufunc.py:151: RuntimeWarning: invalid value encountered in get_gradient_sign
  return self.ufunc(*args, **kwargs)

ed.plot.histogram_turb_estimates(natre)

\(R_ρ, K_T, K_ρ\)#

(
    natre.Rρ.sel(pres=slice(250, None))
    .where(np.abs(natre.Rρ) < 100)
    .mean(["latitude", "longitude"])
    .reset_coords(drop=True)
    .hvplot()
).opts(show_grid=True)

ds = xr.Dataset(
    {"KrhoKt": (natre.Krho / natre.Kt), "Rρ": natre.Rρ.where(np.abs(natre.Rρ) < 5)}
).sel(pres=slice(250, None))
ds.plot.scatter(y="KrhoKt", x="Rρ")
plt.yscale("log")

import colorcet
import datashader as dshdr
from holoviews.operation.datashader import datashade

df = ds.to_dataframe()
# df["KrhoKt"] = np.log10(df.KrhoKt)
pts = hv.Points(df, ["KrhoKt", "Rρ"])
datashade(pts, cmap=colorcet.fire).opts(width=600, height=600)

TS plot#

hdl, axes = dcpy.oceans.TSplot(
    natre.salt[::2],
    natre.theta[::2],
    Pref=1000,
    # rho_levels=bins,
    Sbins=np.arange(35, 35.5, 0.025),
    Tbins=np.arange(4, 10, 0.3),
    hexbin=False,
    plot_kwargs={"alpha": 0.01},
)

(
    natre["chi"].rolling(depth=100, center=True, min_periods=1).mean().compute()
).plot.line(
    y="depth",
    col="longitude",
    row="latitude",
    yincrease=False,
    xlim=[1e-12, 1e-7],
    ylim=(2000, 0),
    xscale="log",
)

tornado.application - ERROR - Exception in callback functools.partial(<bound method IOLoop._discard_future_result of <zmq.eventloop.ioloop.ZMQIOLoop object at 0x7f70c8187340>>, <Task finished name='Task-481' coro=<Cluster._sync_cluster_info() done, defined at /home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py:104> exception=OSError('Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s')>)
Traceback (most recent call last):
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 284, in connect
    comm = await asyncio.wait_for(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/asyncio/tasks.py", line 498, in wait_for
    raise exceptions.TimeoutError()
asyncio.exceptions.TimeoutError

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 741, in _run_callback
    ret = callback()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 765, in _discard_future_result
    future.result()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py", line 105, in _sync_cluster_info
    await self.scheduler_comm.set_metadata(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 785, in send_recv_from_rpc
    comm = await self.live_comm()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 742, in live_comm
    comm = await connect(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 308, in connect
    raise OSError(
OSError: Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s
tornado.application - ERROR - Exception in callback functools.partial(<bound method IOLoop._discard_future_result of <zmq.eventloop.ioloop.ZMQIOLoop object at 0x7f70c8187340>>, <Task finished name='Task-488' coro=<Cluster._sync_cluster_info() done, defined at /home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py:104> exception=OSError('Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s')>)
Traceback (most recent call last):
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 284, in connect
    comm = await asyncio.wait_for(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/asyncio/tasks.py", line 498, in wait_for
    raise exceptions.TimeoutError()
asyncio.exceptions.TimeoutError

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 741, in _run_callback
    ret = callback()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 765, in _discard_future_result
    future.result()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py", line 105, in _sync_cluster_info
    await self.scheduler_comm.set_metadata(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 785, in send_recv_from_rpc
    comm = await self.live_comm()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 742, in live_comm
    comm = await connect(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 308, in connect
    raise OSError(
OSError: Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s

<xarray.plot.facetgrid.FacetGrid at 0x7f70a01c4e50>

Procedure#

The large-scale average operator ⟨⟩ represents

a horizontal average over the survey lateral scale,
a vertical average over O(100) m, and
a time average over the 18-day survey.

The mean fields are derived by averaging all variables along neutral-density surfaces \(γ_n\)

Mean vertical gradients \(∂_z θ_m\), \(∂_z S_m\), and \(∂_z b_m\), are calculated from O(100)-m linear fits to \(θ_m(z_n)\), \(S_m(z_n)\), and \(b_m(z_n)\), where \(z_n\) is the mean depth of each surface \(γ_n\).

large-scale temperature gradient \(|∇_n θ_m|\) is estimated by fitting a plane to the 100 stations in the 400 km × 400 km large-scale survey grid and rms is computed as the standard deviation of the departures from the plane fit.

In theory, these bins are approx. 100m apart in neutral density

Choose bins#

# bins = (
#    natre.gamma_n.mean(["latitude", "longitude"])
#    .interp(depth=np.arange(150, 2001, 100))
#    .dropna("depth")
#    .data
# )

bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
natre.gamma_n.stack({"latlon": ["latitude", "longitude"]}).drop("latlon").cf.plot(
    hue="latlon",
    y="depth",
    lw=0.5,
    add_legend=False,
    yincrease=False,
)
dcpy.plots.linex(bins)

Computing along latitude and longitude lines#

def compute_single_line(section):
    return ed.sections.bin_average_vertical(
        section.cf.stack({"cast": ("latitude", "longitude")}).drop("cast"),
        "neutral_density",
        bins,
    )


lonlines = natre.groupby("longitude", squeeze=False).map(compute_single_line)
latlines = natre.groupby("latitude", squeeze=False).map(compute_single_line)

latlines.to_netcdf("../datasets/natre-along-lat-lines.nc")
lonlines.to_netcdf("../datasets/natre-along-lon-lines.nc")

Buoyancy Reynolds Number#

There are a bunch of low values

natre["ν"] = dcpy.oceans.visc(natre.salt, natre.temp, natre.pres)
natre["κ"] = dcpy.oceans.tdif(natre.salt, natre.temp, natre.pres)
natre["Reb"] = natre.eps / natre.ν / natre.N2
natre["Reb"].attrs["long_name"] = "$Re_b$"
del natre.Reb.attrs["units"]

natre["χmol"] = 2 * natre.κ * natre.Tz**2

_, bins, _ = np.log10(natre.χmol + natre.chi).plot.hist(
    bins=101, histtype="step", density=True, lw=2, aspect=3, size=3
)
np.log10(natre.chi).plot.hist(bins=bins, histtype="step", density=True, lw=2);

kwargs = dict(bins=np.linspace(-2, 5, 101), density=True, histtype="step", lw=2)

np.log10(natre.Reb.sel(latitude=24.7, method="nearest")).plot.hist(
    **kwargs, aspect=3, size=3
)
np.log10(natre.Reb.sel(latitude=25.5, method="nearest")).plot.hist(**kwargs)

dcpy.plots.linex(np.log10([16, 200]))

/Users/dcherian/work/python/xarray/xarray/core/computation.py:734: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
/Users/dcherian/work/python/xarray/xarray/core/computation.py:734: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)

NATRE \(ε\) vs \(ε_χ\)#

np.log10(
    natre.eps_chi.where(np.abs(natre.Tz) > 5e-3).where(natre.eps_chi < 1e-6)
).plot.hist(bins=101, histtype="step", density=True)
np.log10(natre.eps).plot.hist(bins=101, histtype="step", density=True, yscale="linear");

natre.eps_chi.where(np.abs(natre.Tz) > 5e-3).where(natre.eps_chi < 5e-6).mean(
    ["latitude", "longitude"]
).coarsen(pres=100, boundary="trim").mean().plot()
natre.eps.mean(["latitude", "longitude"]).coarsen(
    pres=100, boundary="trim"
).mean().plot(yscale="log")

[<matplotlib.lines.Line2D at 0x16f00e0a0>]

Compute NATRE microscale budget terms#

In theory, these bins are approx. 100m apart in neutral density

natre = ed.natre.read_natre(load=True).drop("depth")
bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
natre.gamma_n.stack({"latlon": ["latitude", "longitude"]}, create_index=False).cf.plot(
    hue="latlon",
    y="pres",
    lw=0.5,
    add_legend=False,
    yincrease=False,
)
dcpy.plots.linex(bins)

micro.load(scheduler=client).to_netcdf("../datasets/natre-density-binned.nc")

OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

T profiles#

ds = natre.stack(
    {"station": ("latitude", "longitude")}, create_index=False
).assign_coords(station=np.arange(100))

import xfilter

Zname = ds.cf.axes["Z"][0]

CT = ds.CT.interpolate_na(Zname)
sortlo = xfilter.lowpass(CT, coord=Zname, freq=1 / 20)
sortlo100 = xfilter.lowpass(CT, coord=Zname, freq=1 / 100)

(
    sortlo100.hvplot.line(y=Zname, color="k")
    # * unsorted.CT.hvplot.line(y=Zname)
    * sortlo.hvplot.line(y=Zname, color="r")
    # * sort.CT.hvplot.line(y="pressure", ylim=(2000, 0))
).opts(ylim=(1200, 800), xlim=(6, 10), frame_width=500, frame_height=400)

Better \(\widetilde{K_T T_z}\) estimate#

It appears that preserving the sign of \(wT\) is important. After a lot of averaging I do get nearly all positive values but the error bars are big!

Yeah the sign for \(T_m + T_e\) can be negative and opposite to \(T_z^m\) Indeed that is important!
If I sort the \(Θ\) then I get rid of the stable inversions that are important.
~No matter what I do, I basically recover χ :(.~ I recover χ if I don’t pay attention to this fact.
I think this explains why \(\widetilde{K_t T_z}\) with sign preserved is closer to \(Γ ⟨ε⟩/N_m^2 T_z^m\). The cancellations are important! Using absolute values kills the cancellations…
Might be better to stabilize the sign somehow; Logic here is that
- the sign of \(\widetilde{T_z}\) is set by the intermediate scale (30-50m-ish);
- but we want \(\widetilde{w_tθ_t} = \widetilde{χ}/\widetilde{T_z}/2\) over 5-m scales
- So filter at 10m wavelength; choose sign;
- Then determine \(\widetilde{T_z}\) magnitude over 5-m linear fits;
- apply sign of first element in window (this is somewhat arbitrary)

natre = ed.natre.read_natre(load=True)

natre_binned = xr.load_dataset("../datasets/natre-density-binned.nc")

Tzmi = natre_binned.dTdz_m.reset_coords(drop=True).interp(gamma_n=natre.gamma_n)

Tmi = natre_binned.theta.reset_coords(drop=True).interp(gamma_n=natre.gamma_n)

Depth-space#

Here I estimate \(\widetilde{w^tθ^t} = \widetilde{χ} / 2 / \widetilde{T}_z\) over 5m scales.

Using linear fit to \(T\) (or sorted \(T\)) ofr gradient;
simple 5m average of χ

The core problem is that a simple averaging of \(\widetilde{T}_z\) basically yields \(T_z^m\); but it should be bigger for this strategy to work.

(natre.Tz - Tzmi).plot.line(row="latitude", col="longitude")

<xarray.plot.facetgrid.FacetGrid at 0x17cf48d60>

import hvplot.xarray
import xfilter

loc = dict(latitude=3, longitude=8)

T = natre.CT.isel(loc).interpolate_na("pres").dropna("pres")
# Tm = xfilter.lowpass(T, coord="pres", freq=1 / 200, order=1)
Te = xfilter.lowpass(T, coord="pres", freq=1 / 20, order=1)

# (T.plot())
# (Te.plot())
# (Tmi.isel(loc).plot())
# xfilter.lowpass(
#    (T - Tmi.isel(loc)).interpolate_na("pres"), coord="pres", freq=1 / 20
# ).plot()
# (Tm.plot())
# ((T - Te).plot(ax=plt.gca().twinx()))

# (T.hvplot.line(x="pres") * Te.hvplot.line(x="pres"))
(T.differentiate("pres").hvplot.line()) * Te.differentiate("pres").hvplot.line()

clean.dTdz.isel(latitude=2, longitude=6).plot()
clean.dTdz.isel(latitude=9, longitude=8).plot()
natre_binned.dTdz_m.plot(x="pres", color="k")

[<matplotlib.lines.Line2D at 0x17c902df0>]

(clean.dTdz).coarsen(
    {"latitude": 10, "longitude": 10, "pres_5.1": 20}, boundary="trim"
).mean().plot(yscale="log")
(natre_binned.dTdz_m).plot(x="pres")

[<matplotlib.lines.Line2D at 0x17757ff40>]

clean = ed.sections.estimate_microscale_stirring_depth_space(
    natre, filter_len=10, segment_len=3
)

plt.figure()
clean["Tz~"].coarsen(pres=200, latitude=10, longitude=10, boundary="trim").mean().plot()
natre_binned.dTdz_m.plot(x="pres", color="k", yscale="log")

plt.figure()
natre.KtTz.where(np.abs(natre.Tz) > 1e-5).mean(["latitude", "longitude"]).plot(
    yscale="log", lw=0.5, label="simple"
)
(
    (
        natre.KtTz.where(np.abs(natre.Tz) > 1e-5)
        .coarsen(pres=200, latitude=10, longitude=10, boundary="trim")
        .mean()
    ).plot(yscale="log", label="simple coarsened")
)

(
    clean["KtTz~"]
    # .where(np.abs(clean["KtTz~"]) < 5e-6)
    .cf.coarsen({"Z": 200, "latitude": 10, "longitude": 10}, boundary="trim").mean()
    # .mean(["latitude", "longitude"])
    .plot(label="cleaner")
)

(natre_binned.KρTz2 / natre_binned.dTdz_m).plot(label="$K_ρT_z^m$", color="k", x="pres")
(natre_binned.chib2 / natre_binned.dTdz_m).plot(
    label="$χ/2/T_z^m$",
    color="k",
    x="pres",
    ylim=(1e-8, 1e-6),
)
# (micro.wTTz / micro.dTdz_m).plot(label="$K_tT_z$", color="cyan", x="pres")
plt.legend()

7 20
50243  values < 0

<matplotlib.legend.Legend at 0x177280df0>

kwargs = dict(pres=100, latitude=10, longitude=10, boundary="trim")

clean.KtTz.mean(["latitude", "longitude"]).coarsen(
    {"pres_5.1": 20}, boundary="trim"
).mean().plot()

(
    ((natre.KtTz).where(np.abs(natre.Tz) > 1e-3).coarsen(**kwargs).mean()).plot(
        yscale="log"
    )
)

((natre.KrhoTz).coarsen(**kwargs).mean().plot(yscale="log", color="k"))
# np.abs(natre.chi / 2).mean(["latitude", "longitude"]).coarsen(**kwargs).mean().plot(
#    yscale="log"
# )

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [316], in <cell line: 3>()
      1 kwargs = dict(pres=100, latitude=10, longitude=10, boundary="trim")
----> 3 clean.KtTz.mean(["latitude", "longitude"]).coarsen(
      4     {"pres_5.1": 20}, boundary="trim"
      5 ).mean().plot()
      7 (
      8     (((natre.KtTz).where(np.abs(natre.Tz) > 1e-3).coarsen(**kwargs).mean())).plot(
      9         yscale="log"
     10     )
     11 )
     13 ((natre.KrhoTz).coarsen(**kwargs).mean().plot(yscale="log", color="k"))

NameError: name 'clean' is not defined

natre.eps.coarsen(**kwargs).mean().plot(yscale="log")
natre.eps_chi.where(np.abs(natre.Tz) > 1e-3).coarsen(**kwargs).mean().plot(yscale="log")

[<matplotlib.lines.Line2D at 0x16469f1f0>]

Interestingly excluding small gradients excludes small χ, so the \(\widetilde{K_TT_z}\) estimate increases with more filtering!

kwargs = dict(bins=100, density=True, histtype="step")
np.log10(chi / 2 / dTdz).plot.hist(**kwargs)
np.log10(chi / 2 / dTdz.where(dTdz > 5e-3)).plot.hist(**kwargs);

T-space#

This is not as easy.

Choosing T bins is a real pain. Because we are covering a giant region; the spread in cast-mean T is large. So picking bins is a little challening.
By averaging over all casts I end up with \(⟨χ⟩\), and without accounting for sign of \(Δz\) (separation between isotherms), I end up with \(T_z^m\) in the denomoinator.

natre = ed.natre.read_natre(load=True)

natre = natre.update(clean)

from xarray.tests import raise_if_dask_computes

bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
with raise_if_dask_computes():
    micro = ed.sections.bin_average_vertical(
        natre.cf.stack(
            {"cast": ("latitude", "longitude")}, create_index=False
        ).assign_coords({"cast": ("cast", np.arange(100), {"cf_role": "profile_id"})}),
        "neutral_density",
        bins,
        blocksize=10,
    )
# micro.pres.attrs["bounds"] = "pres_err"
micro

<xarray.Dataset>
Dimensions:         (bound: 2, gamma_n: 18, cast: 100)
Coordinates:
  * bound           (bound) <U5 'lower' 'upper'
  * gamma_n         (gamma_n) float64 26.56 26.78 26.95 ... 27.93 27.95 27.96
    num_obs         (gamma_n) int64 19358 19463 20127 ... 20191 11530 22215
    pres            (gamma_n) float64 201.4 300.6 401.1 ... 1.795e+03 1.877e+03
    gamma_n_bounds  (bound, gamma_n) float64 26.43 26.69 26.87 ... 27.95 27.97
    latitude        (cast) float64 27.5 27.5 27.5 27.5 ... 23.9 23.9 23.9 23.9
    longitude       (cast) float64 -30.7 -30.3 -29.8 -29.4 ... -27.6 -27.2 -26.8
    time            (cast) datetime64[ns] 1992-03-28T15:28:59.999996928 ... 1...
  * cast            (cast) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
Data variables: (12/54)
    chi             (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    eps             (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    KtTz            (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    KtTz~           (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    hm              (gamma_n) float64 98.77 98.73 101.8 ... 101.3 58.27 114.8
    wT              (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
    ...              ...
    δKtTz~Tz        (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    δwTTz           (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
    δresidual       (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    δpres           (gamma_n) float64 49.38 49.37 50.92 ... 50.63 29.14 57.4
    chib2           (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    chib2_err       (gamma_n, bound) float64 dask.array<chunksize=(1, 2), meta=np.ndarray>
Attributes: (12/14)
    Conventions:           CF-1.6
    netcdf_version:        4
    project:               North Atlantic Tracer Release Experiment (NATRE)
    expocode:              32OC250_4
    cast_number:           3.0
    title:                 Microstructure profiler data from the ship Oceanus...
    ...                    ...
    longitude:             -30.723333333333333
    chief_scientist:       Raymond W. Schmitt
    data_originator:       Polzin
    institution:           WHOI
    data_assembly_center:  CCHDO
    commit:                eddydiff: 3c5e37fac5bc4c53f9567addc4c058ef415a63b4...

Using fancy sign detection#

def notnull(data):
    flat = data.data.ravel()
    return flat[~np.isnan(flat)]

mask = np.abs(clean["KtTz~"]) > 1e-5
clean["chi~"].where(mask).pipe(notnull)

array([2.29729636e-05, 4.99936064e-06, 8.38180596e-08, 4.42592055e-08,
       2.45233000e-08, 6.14768901e-08, 5.21607218e-07, 2.45605935e-07,
       9.83617336e-08, 2.11507639e-07, 1.67454668e-08])

clean["Tz~"].where(np.abs(clean["KtTz~"]) > 1e-5).pipe(notnull)

array([ 0.14074649,  0.13337417,  0.00319876,  0.00175999,  0.00105803,
        0.00300432,  0.00974087,  0.00935505, -0.00478368,  0.00899659,
        0.00078514])

clean["Tz~"].isel(cast=16).plot(marker="x")

[<matplotlib.lines.Line2D at 0x191bccb80>]

clean["KtTz~"].isel(cast=26).ffill("pres").plot(yscale="linear")

[<matplotlib.lines.Line2D at 0x193737850>]

ed.plot.debug_section_estimate(micro)

dcpy.plots.fill_between_bounds(micro, "chib2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KρTz2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KtTz~Tz", y="pres")
# dcpy.plots.fill_between_bounds(micro, "KtTzTz", y="pres")
# dcpy.plots.fill_between_bounds(micro, "wTTz", y="pres")
plt.xlim([1e-11, 1e-8])
plt.xscale("log")

Single bin#

from xarray.tests import raise_if_dask_computes

bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
with raise_if_dask_computes():
    grouped = ed.sections.bin_average_vertical(
        natre.cf.stack(
            {"cast": ("latitude", "longitude")}, create_index=False
        ).assign_coords({"cast": ("cast", np.arange(100), {"cf_role": "profile_id"})}),
        "neutral_density",
        bins[10:],
        blocksize=10,
        return_group=True,
    )
# micro.pres.attrs["bounds"] = "pres_err"
for label, group in grouped:
    break

Tbins, sortT = ed.sections.estimate_microscale_stirring(
    group.unstack(), dz=7, debug=True
)
dT = np.diff(Tbins).mean()
sort = dcpy.oceans.thorpesort(sortT, by="gamma_n", core_dim="pres")
sort.CT.plot.line(hue="cast", add_legend=False);

sort.CT.isel(cast=[20, 30, 60]).cf.dropna("Z").plot(hue="cast")
dcpy.plots.liney(Tbins)

CT = profile.CT.data

profiles = natre.copy(deep=True)
profiles["CT"] = profiles.CT.interpolate_na("pres")
profiles["Tzfilt"] = -1 * xfilter.lowpass(
    profiles.CT, coord="pres", freq=1 / 10
).differentiate("pres").bfill("pres")

stacked = profiles.stack({"cast": ("latitude", "longitude")}, create_index=False)

profile = stacked.interpolate_na("pres")  # .sel(pres=slice(1120, 1250))


# get_gradient_sign(
#    profile.pres.data, profile.CT.data, profile.Tzfilt.data, 5, 10, a.data
# )
# profile.CT.plot()
Tzsign.plot(ax=plt.gca().twinx())

<matplotlib.collections.QuadMesh at 0x185d14790>

def clean_idx(idx, CT, ΔP, window):

    P = CT["pres"].data
    CT = CT.data

    def _clean(idx, P, ΔP, window):
        newidx = [idx[0]]
        for left, right in zip(idx[:-1], idx[1:]):
            Pl, Pr = P[left], P[right]
            # print(np.abs(Pr - Pl))
            if np.abs(Pr - Pl) > ΔP:
                if ((right - left) % window) < window // 2:
                    right += 1
                newidx.append(right)

        if newidx[-1] != idx[-1]:
            # Last index was dropped
            # instead to preserve all data;
            # we drop the previous one and re-add the last one
            newidx.pop()
            newidx.append(idx[-1])
        print(idx, newidx)
        return newidx

    idx = _clean(idx, P, ΔP, window=window)
    idx = _clean(idx, CT, 1e-4 * ΔP, window=1e10)
    return idx


def smooth_gradient(profile, chi, ΔP, window, debug=False):

    dp = 0.5
    smooth = xfilter.lowpass(
        profile, coord="pres", freq=1 / int(ΔP // dp), num_discard=5
    )
    Tzsmooth = smooth.differentiate("pres")
    # array = profile.data
    if debug:
        f, ax = plt.subplots(4, 1, sharex=True, constrained_layout=True)
        profile.plot(ax=ax[0])
        smooth.plot(ax=ax[0])

    (idx,) = np.nonzero(np.abs(np.sign(Tzsmooth).diff("pres").data) > 0)
    idx = np.concatenate([[0], idx, [profile.sizes["pres"] - 1]])
    if debug:
        dcpy.plots.linex(Tzsmooth["pres"][idx], color="k", lw=0.5, ax=ax[0])

    idx = clean_idx(idx, profile, ΔP=ΔP, window=window)
    if debug:
        dcpy.plots.linex(Tzsmooth["pres"][idx], color="b", ax=ax[0])

    Tzs = []
    for subprof in np.split(profile, idx[1:-1]):
        # -1 for temperature slope
        sign = np.sign(subprof[0] - subprof[-1]).data
        coarse = subprof.coarsen(pres=window, boundary="pad").construct(
            pres=("p", "window")
        )
        sorty = dcpy.oceans.thorpesort(
            coarse, coarse, core_dim="window", ascending=False
        )
        sorty = sorty.where(sorty.count("window") > 2)
        Tzsubprof = sign * np.abs(
            sorty.polyfit("window", deg=1)
            .sel(degree=1, drop=True)
            .polyfit_coefficients.data
        )
        Tzs.append(np.repeat(Tzsubprof, window)[..., : subprof.shape[-1]])

    dTdz = profile.copy(data=np.concatenate(Tzs)).ffill("pres")
    if debug:
        dTdz.plot(ax=ax[1])
        (-1 * Tzsmooth).plot(ax=ax[1])

        (chi).plot(yscale="log", ylim=(1e-12, None), ax=ax[2])
        ax2 = ax[2].twinx()
        (dTdz).plot(ax=ax2, color="k")
        (dTdz.where(np.abs(dTdz) < 2e-4)).plot(color="b", marker="o", ls="none", ax=ax2)
        dcpy.plots.set_axes_color(ax2, "k")

        (1025 * 4000 * -chi / dTdz).plot(ax=ax[3], yscale="symlog")
        (1025 * 4000 * -chi / -Tzsmooth).plot(ax=ax[3])
        (1025 * 4000 * -chi / dTdz.where(np.abs(dTdz) < 2e-4)).plot(
            color="b", marker="o", ls="none", ax=ax[3]
        )

        dcpy.plots.clean_axes(np.atleast_2d(ax).T)
        f.set_size_inches((10, 7))

    return dTdz
    # dcpy.plots.liney(Tbins)


# NATRE 36, 56, 60 are good to test with
profile = sort.isel(cast=36).cf.dropna("Z")

dTdz = smooth_gradient(profile.CT, profile.chi, ΔP=5, window=20, debug=True);

[  0   6  11  37  59  76 104 119 124 135 148 153 187] [0, 38, 60, 76, 105, 119, 135, 148, 187]
[0, 38, 60, 76, 105, 119, 135, 148, 187] [0, 39, 61, 77, 106, 120, 136, 149, 187]

Tz_m = sort.CT.polyfit("pres", deg=1).sel(degree=1, drop=True).polyfit_coefficients
Tz_m.plot()

[<matplotlib.lines.Line2D at 0x180d0ceb0>]

Tz = -1 * xfilter.lowpass(sort.CT, coord="pres", freq=1 / 10).differentiate("pres")
KtTz = sort.chi / Tz
KtTz = KtTz.where(np.abs(KtTz) < 300 / 1025 / 4000)

(KtTz.mean("pres") * Tz_m).plot()
(sort.chi.mean("pres") / 2).plot()

[<matplotlib.lines.Line2D at 0x1814e42b0>]

Try ε#

natre.eps_chi.where(np.abs(natre.Tz) > 1e-10).mean(["latitude", "longitude"]).plot(
    yscale="log"
)
natre.eps_chi.where(np.abs(natre.Tz) > 1e-2).mean(["latitude", "longitude"]).plot(
    yscale="log"
)
natre.eps.mean(["latitude", "longitude"]).plot(yscale="log")

[<matplotlib.lines.Line2D at 0x170f51220>]

Sensitivity of χpod estimate to mean K estimate#

Various estimates of the variance produced by turbulent stirring of mean vertical gradient

Ferrai & Polzin (2005): \(Γ ⟨ε⟩/N_m^2 ∂_zθ_m^2\)
χpod estimate: \(⟨K_T θ_z⟩ ∂_zθ_m\)
(1) with \(ε_χ\) : \(Γ ⟨ε_χ⟩/N_m^2 ∂_zθ_m^2\)
(2) with \(K_ρ\) instead of \(K_T\): \(⟨K_ρ θ_z⟩ ∂_zθ_m^2\)

(chidens.Krho_m * chidens.dTdz_m**2).cf.plot.step(
    color="k", lw=2, label="$Γ ⟨ε⟩/N_m^2  ∂_zθ_m^2$"
)

(chidens.KtTz * chidens.dTdz_m).cf.plot.step(
    color="k", marker="x", ls="none", label="$⟨K_T θ_z⟩ ∂_zθ_m$"
)

(0.2 * chidens.eps_chi / chidens.N2_m * chidens.dTdz_m**2).cf.plot.step(
    color="C1", lw=2, label="$Γ ⟨ε_χ⟩/N_m^2  ∂_zθ_m^2$"
)

(chidens.KrhoTz * chidens.dTdz_m).cf.plot.step(
    lw=2, label="$⟨K_ρ θ_z⟩  ∂_zθ_m$", color="C2"
)

plt.grid(True, which="both", lw=0.5)

plt.legend()
plt.xscale("log")
plt.xlabel("Variance production or dissipation [°C²/s]")
plt.gcf().set_size_inches((4, 5))

f, ax = plt.subplots(1, 2, sharey=True, sharex=True)
micro.Krho.cf.plot.step(y="pres", xscale="log", ax=ax[0])
micro.Kt.cf.plot.step(y="pres", ax=ax[0])
micro.Krho_m.cf.plot.step(y="pres", ax=ax[1])
micro.Kt_m.cf.plot.step(y="pres", ax=ax[1])
ax[0].legend(["$⟨K_ρ⟩$", "$⟨K_T⟩$"])
ax[1].legend(["$K_ρ^m = Γ ⟨ε⟩/N_m²$", "$K_T^m = ⟨χ⟩/2θ_m^2$"])
[aa.grid(True, which="both", lw=0.5) for aa in ax]
plt.gcf().set_size_inches((6, 5))

Eddy Diffusivity

NATRE microstructure dataset

Contents